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Abstract 

We have developed and implemented a numerical evolution scheme for a class of stochastic prob- 
lems in which the temporal evolution occurs on widely-separated time scales, and for which the 
slow evolution can be described in terms of a small number of moments of an underlying probability 
distribution. We demonstrate this method via a numerical simulation of chemotaxis in a popula- 
tion of motile, independent bacteria swimming in a prescribed gradient of a chemoattractant. The 
microscopic stochastic model, which is simulated using a Monte Carlo method, uses a simplified de- 
terministic model for excitation/adaptation in signal transduction, coupled to a realistic, stochastic 
description of the flagellar motor. We show that projective time integration of "coarse" variables 
can be carried out on time scales long compared to that of the microscopic dynamics. Our coarse 
description is based on the spatial cell density distribution. Thus we are assuming that the system 
"closes" on this variable so that it can be described on long time scales solely by the spatial cell 
density. Computationally the variables are the components of the density distribution expressed 
in terms of a few basis functions, given by the singular vectors of the spatial density distribution 
obtained from a sample Monte Carlo time evolution of the system. We present numerical results 
and analysis of errors in support of the efficacy of this time-integration scheme. 
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I. INTRODUCTION 



A recurring bottleneck in computational modeling of physical processes is the existence 
of multiple space and/or time scales. Important examples include patterns in fluids, defect 
dynamics in solids, and molecular dynamics of macromolecules. For example, in the latter 
case, the evolution of the physical configuration from initial to final state typically occurs 
on the time scale of milliseconds, while the interactions between constituent components 
must be resolved on the picosecond time scale. Furthermore, it is increasingly important 
in simulating biological systems to have a hybrid computational framework for interact- 
ing stochastic and deterministic processes, which occur on different time scales. Multiscale 
computational methods combine processing on fine space/time scales according to the gov- 
erning^microscopic description with macroscopic changes on a coarse grid (see for example, 

mm 

Here we present a computational scheme for "coarse projective" time integration of the 
macroscropic dynamics of stochastic processes that involve multiple, widely-separated time 
scales. To illustrate the method we apply it to a microscopic model that generates the 
biased random walk describing the macroscopic motion of a bacterium such as E. coli in 
an attractant gradient. Details of this Monte Carlo model of bacterial chemotaxis that 
integrates signal transduction with the motor response are given in Section ITTTl The focus 
of our analysis in this paper is the computational scheme, and as such, our results address 
its plausibility rather than the phenomenology of chemotaxis. 

In Section|n]we discuss previous work on coarse time integration of microscopic dynamics 
and the new features of the method developed here, and we motivate its application to 
chemotaxis. The Monte Carlo description is presented in Section UTTl We outline the coarse 
integration scheme in Section IIVI and present numerical results and analysis of errors in 
Section El We conclude with future extensions and applications of this work. 

II. PREVIOUS WORK 

The traditional approach to studying long-term dynamics of multiscale processes involves 
(a) the derivation of a "coarse-grained" set of evolution equations, followed by, (b) the ana- 
lytical and/or computational study of these reduced equations using established, continuum 
numerical analysis tools. Recently, a so-called "equation-free" approach to the study of the 
coarse-grained behavior of such problems has been proposed which circumvents the first 
step (USE])- This computational approach is based on the "coarse," or macroscopic time- 
stepper, a map from the coarse variables at time t = to those at time t = T, where T is 
typically much larger than characteristic microscopic time scales in the system. This map 
is not obtained directly from the macroscopic evolution equations, which we may not know, 
but rather through short time evolution intervals of appropriately initialized microscopic 
simulations. The initial macroscopic variables are lifted to microscopic variables to initialize 
the microscopic simulation. At the completion of a burst of microscopic simulation, the mi- 
croscopic variables are restricted back to macroscopic variables, providing an approximation 
to the macroscopic time step. This provides a chord of the macroscopic solution, which is 
an approximation to the time derivative of the macroscopic solution. This value can then 
be used in any conventional continuum numerical method for the macroscopic equations. 
This approach has been apphed in several microscopic contexts, and the results appear to 
be promising 0, H, H, [l^ . 
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In this work we apply the coarse timestepper in a projective integration study of a 
spatially-distributed kinetic Monte Carlo simulation of a biased random walk. When, as 
in this case, the microscopic equations are stochastic, an effective algorithm must reduce 
the variance inherent in individual realizations of the stochastic process to a level that can 
be tolerated by the continuum numerical algorithm applied to the coarse system. This can 
be done by lifting to multiple copies, or, as is done here, by estimating the derivative from 
a least-squares fit to a large number of microscopic time steps. The projective integration 
method applied here uses a derivative estimate for the projective step, so that the stochastic 
noise is amplified by the inverse of the effective step length used in the derivative estimate. 
Hence, variance reduction is very important. Augmenting the number of copies of the sim- 
ulation (which, for our noninteracting particle model corresponds also to a simulation with 
a larger number of cells) is the most direct approach to variance reduction; other variance 
reduction schemes are discussed in (ll . 1^ . 

In the context of bacterial chemotaxis, the computational scheme developed here can 
be viewed as a direct method to compute the macroscopic evolution of the cell density in 
space and time without actually deriving these equations. An alternate approach begins 
with the transport equation for the velocity jump process, in which discontinuous changes 
in the direction (or speed) of an individual cell are generated by a Poisson process. It can 
be shown rigoro usly that this reduces to the chemotaxis equation under suitable scaling of 
space and time jl3| 

^ = V-(DV/io-/iox(5)V5), (1) 

In the above, fio{x,t) represents the density of particles at spatial position x at time t, D 
is the diffusion constant, S is the concentration of the chemotactic attractant /repellent, 
and x{S) is the chemotactic sensitivity. However, there is as yet no analytical procedure 
available for determining the diffusion constant and chemotactic sensitivity when dependence 
on internal state variables determining the cell's response to the external signal is explicitly 
included the transport equation. We show that a macroscopic time-stepper can be used here 
in the time evolution of the spatial density, even though the macroscopic evolution equations 
are not known. 



III. BACTERIAL CHEMOTAXIS 
A. Signal transduction 

The ability to sense and respond to environmental cues is necessary for the survival of 
most organisms. Escherichia Coh {E. coli) is a common and well-known single cell organ- 
ism, with roughly 4000 genes, whose chemotaxis network has emerged as a prototype for 
understanding signal transduction networks in general 3|. Its genome is known, the crystal 



structures of many proteins have been obtained, and a large number of mutant strains exist, 
allowing detailed behavioral studies. 

For each cell, chemotactic behavior begins when attractant or repellent molecules bind 
to membrane receptors, triggering a cascade of chemical reactions inside the cell that culmi- 
nates in the production of the phosphorylated form of a response regulator protein (CheY-P) 
which controls the direction of rotation of the flagellar motor. The series of reactions that 
converts the extracellular signal (attractrant/repellent) into cellular response is referred to 
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as the signal transduction pathway. Flagella possess an inherent chirahty, such that counter- 
clockwise (CCW) rotation results in bundling-up of the six to eight flagella per cell, allowing 
them to act as a single propeller and leading to smooth swimming motion of the bacterium. 
When the flagella rotate clockwise (CW), the bundle flies apart and the bacterium tumbles. 
As conditions become increasingly favorable due to increase in chemoattractant concentra- 
tion, a bacterium extends its run-length; otherwise, it tumbles, and the subsequent direction 
of motion is randomly chosen, allowing a more favorable direction to be discovered. The 
resulting motion is a biased random walk toward favorable conditions and away from less 
favorable ones. 

The sequence of biochemical reactions that take place inside a bacterium, starting with the 
binding of an attractant or repellent molecule to receptors on the cell surface and leading to 
the change in concentration of the response regulator species, CheY-P, has been extensively 
studied. Both deterministic models of these biochemical pathways, based on the law of mass 
action and Michaelis-Menten kinetics jisl . 3], as well as fully stochastic models exist jl7l |. 



Measurement of concentration change is achieved through a temporal rather than a spatial 
comparison: fast sampling of the present external concentration is compared with memory 
of that concentration some time ago. Order-of-magnitude analyses for why measurement of 
concentration ch ang es as spatial gradients across the cell length is not physically feasible 
have been given [ISf. 

Memory in the signal transduction network is achieved through the existence of fast and 
slow reaction time scales. The fast reactions are receptor-ligand binding and phosphorylation 
kinetics; the slow reactions are methylation and demethylation. Ligand binding reduces the 
autophosphorylation rate of the corresponding membrane-bound receptor, in turn decreasing 
the rate of transfer of phosphoryl groups to CheY, and resulting in a (fast) drop in [CheY- 
P]. Addition of methyl groups to ligand-bound receptors restores the autophosphorylation 
rate, resulting in (slow) increase in [CheY-P]. The rate of demethylation becomes significant 
once a high methylation level is achieved. Hence, [CheY-P] reflects the balance between the 
fraction of ligand-bound receptors and methylated receptors. Perfect adaptation refers to the 
return of [CheY-P] to the same steady state level, regardless of the constant concentration 
level of the external stimulus. This value falls within the fixed operational range of the 
motor response to CheY-P. 

A minimal model representing the signal transduction process includes fast excitation 
and slow adaptation to an external stimulus and is given in [l9| : 

dui _ {f{S) - U2) - ui 

dt ~ r, ' ^> 

dU2 _ f{S) - U2 , s 

dt - ' 

where Te and Ta are the excitation and adaptation times, respectively, with Te -C r^. We 
identify Ui with the deviation of [CheY-P] from its steady state value, and U2 as the number 
of methylated receptors per unit volume. f{S) = f{S(t)) is a function of the external 
stimulus; for example, it is the number of bound receptors per unit volume in the presence 
of an external signal concentration, S: 

f = NTS/{KL + S). (4) 

Nj- is the total number of receptors per unit volume, and is the dissociation constant 
for the ligand-binding process. 
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The formal solution to Eqs. (121)^© can be easily obtained for each cell as 



u,{t) = ^ t f{S{t'))//^-dt', (5) 




u,{t) = '-^i f f{S{t'))e''/^^dt' r e-t'{i/ra-i/r.) C f^S{t"))e'"'^-dt"dA{<6) 



However, because S{t) = S{x(t),t), where x{t) is a random variable that represents the cell's 
position at time t in a given external concentration field, the integration takes place along the 
cell trajectory, which is a biased random walk: x{t') for t' < t biases the probability of x{t). 
Hence, S(t) is a stochastic variable, and the integration must be carried out computationally. 

B. Motor and cell response 

We adopt a stochastic approach to modeling the response of the flagellar motor to the 
regulator species, CheY-P. The motor is assumed to be in one of two states, CW or CCW, 



corresponding to clockwise or counterclockwise rotation j20|, |21|. The rates of transition 
between these states are functions of [CheY-P]: 



fc4 



CCWi CW, t = l,...,Nj, (7) 



k 



where Nf is the number of fiagella per cell. For any concentration of CheY-P, the motor has 
a nonzero probability of being in either state. The probability of waiting at least a time T 
for a transition from CCW to CW, $+, or a transition from CW to CCW, is given by: 

If the motor is probed at a time T = At, such that k±At <^ 1, then the probability that a 
change in direction will occur in (T, T + At) can be approximated as 

1 - $± ^ A;±(0)At. (9) 

Systematic studies have been undertaken that show simulation results are qualitatively in- 
dependent of the choice of At once the above restrictions are met. 

In this work, we use equilibrium transition rates, k±, consistent with recent experimental 
measurements [2^ of the single flagellum clockwise bias, feciy, and the reversal frequency, 
w, as functions of [CheY-P]. The rates are related to these measured quantities according 

to ET 



Jew 



(10) 



k.+k 



2k. k^ , , 

where the reversal frequency is the geometric mean of the transition rates. We take hew 
to be described by a Hill function, hew = Vv /{Km^ + yj"), where = [CheY - P], with 
experimentally measured values for the Hill coefficient and dissociation constant, equal to 
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h = 10.3 and Km = 3.1 /iM. In these experiments, it was noted that the data were consistent 
with w = {2fiMs~^) dbcw/dyp. Based on these results, we use 

- ^y^'" (12) 

^■^1 ^f-\. (13) 

In a constant chemoattractant field, Hi, HI show that the "Voting Hypothesis" is suc- 
cessful in producing the correct running bias of the cell, Bccw, from the individual fiagellar 
bias, bccw- 

Bccw = J2[n ^ccw (1 - hccwff-' , (14) 

where u is the minimum number of flagella required to be in the CCW state for the cell 
to run. Although cooperativity among the fiagella is not yet fully understood, we similarly 
adopt the "majority rules" algorithm to determine whether each cell will run or tumble in 
the presence of a chemoattractant gradient: 



V = 

1=1 



J2s^^K (15) 



where s^^^ = 0, 1 is the state of the i^^ fiagellum, corresponding to CW and CCW directions 
of rotation, respectively. Based on the experimentally measured gain of a single fiagellum 
|2^ . a steady state value of ijp = 2.95 fiM yields bccw ~ 0.64 to within experimental error. 
The resulting running bias for the cell, according to Eq. HH is Bccw ~ 0.86 for Nj = 6 
and z/ = 3, in agreement with experimental results. Hence, if three or more fiagella are 
determined to be in the CCW state, the cell in our simulation runs; otherwise, it tumbles. 



C. Monte Carlo scheme 

For each cell, the Monte Carlo scheme can be outlined as follows. At t„: 

(a) Each fiagellum has a state Sn, where = if CW, and s„ = 1 if CCW. To determine 
s„+i, for each fiagellum, we draw a uniformly distributed random number, ( G [0, 1]. 

(b) We determine an approximation to the probability for changing direction of rotation, 
given by Eq. El If the values of k± are such that this approximation is not valid. At is 
reduced until it is. This probability is a function of [CheY-P], which depends on the 
trajectory that the cell has taken along the external chemoattractant gradient. 

— If the fiagellum is in the CW state (s^ = 0) and C, < k_ At, it continues in the 
CW state (s„+i = 0); else, it switches to the CCW state (s„+i = 1). 

— If the fiagellum is in the CCW state (s„ = 1) and ( < k^ At, it continues in the 
CCW state {sn+i = 1); else, it switches to the CW state {sn+i = 0). 
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(c) If three or more flagella are now in the CCW state, Un+i > 3 in Eq. and the cell 
runs; else, it tumbles. If the cell is determined to run and z/„ < 3, the direction of 
motion is chosen to be left or right with equal probability. Otherwise, it continues to 
run in the same direction. 

(d) The position of the j^^ cell, Xn{j), is accordingly updated to Xn+i (j) (using the accepted 
time step At), with the cell speed, Vceii, assumed to be constant and independent of 
position. 

(e) The signal transduction variables (mi,M2) are integrated in time, using the forward 
Euler scheme. Their time history is a function of each cell's trajectory. 



D. Combining signal transduction and motor response 

Although we use a toy excitation-adaptation model to describe signal transduction, it is 
coupled to an experimentally realistic description of the flagellar switch. To do so requires 
introducing (i) a shift in the steady state value of ui, and (ii) an amplification or gain factor, 
go, for the signal transduction network, required so that the output of this network spans 
the dynamic range of the motor: 

yp = Vp- 9oUi (16) 
The issue of gain in chemotaxis continues to receive much attention 2^. Given the level 



of simplicity of our signal transduction model, we do not imply this treatment of network 
gain to be physically realistic. An obvious drawback is that, like other parameters in the 
toy model (excitation and adaptation times), it is static, whereas in the actual biological 
system, these parameters depend on the external signal. 

In Table HJ we report the numerical values for the MC model parameters that operationally 
interface the toy signal transduction model with the realistic model of the flagellar switch. 
Hence, although the choice of toy model parameters was "physically motivated" , we do not 
attempt to make direct correspondence with experimental values. We caution the reader 
against making quantitative connection between our numerical results at these values and 
experimental results. 



TABLE I: Numerical values for parameters used in the Monte Carlo evolution. 



Parameter 


Value 


Ta 


100 s 


Te 


0.1 s 


90 


5 


Nt 


15 ijM 


Kl 


1 //M 


Vp 


2.95 pM 


Km 


3.1 


Vcell 


20 /um s^^ 
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IV. COARSE INTEGRATION 



For the present problem of simulating the chemotactic response of a population of inde- 
pendent bacteria, the microscopic phase space is a directproduct of the phase space for each 
cell, consisting of its 

• signal transduction variables: u={ui,U2), 

• flagellar state: s, where s^*^ = 0, 1 for i = 1 . . . Nf, 

• position: x, 

• direction of motion, d = R,L,T, corresponding to running right, running left and 
tumbling, respectively. 

Rather than evolving all microscopic degrees of freedom, a coarse integration scheme iden- 
tifies suitable reduced variable(s) to be integrated. A clearly relevant reduced variable in 
population studies is the spatial density distribution of cell positions, fio{x, t), obtained from 
the set of individual cell positions, xj, which are a subset of the full microscopic phase space. 
However, it is possible that the unknown equations of the coarse description use additional 
variables, for example, the densities of the right moving, left moving, and tumbling cells. 
Pi where i = R,L, or T. In this case, fio{x,t) = J2iPii^y'^)- Indeed, an optimal reduced 
representation balances accuracy and efficiency of modeling the physical process. Here, we 
propose to retain only the spatial density distribution /io(x, t) as the relevant macroscopic 
variable in coarse modeling of chemotaxis since, as we will demonstrate numerically, the 
density of cells in each state rapidly approaches a functional of the total density. 

A systematic approach to testing the adequacy of a particular coarse description is dis- 
cussed in [9!]. It is based on locating the same fixed point at different levels of coarse descrip- 
tion, and examining the slow eigenvalues and corresponding eigenvectors of the linearization 
of the coarse timestepper at these fixed points (see also the discussion in Q). 

The separation of characteristic time scales describing the underlying macroscopic and 
microscopic dynamics - for example, the (long) time scale on which the spatial density dis- 
tribution of a population of bacteria moving in a chemoattractant profile changes versus the 
(short) mean runtime of a single cell - allows taking time steps in evolving the macroscopic 
state that are "long" in comparison with the "fast" microscopic time scales in the prob- 
lem, resulting in a computationally - efficient time evolution of the macroscopic state. We 
demonstrate this assumption of the separation of time scales explicitly in Section Figure 
^ shows a cartoon sketch of the coarse integration (CI) procedure, where the solid trajec- 
tories denote the restriction of the full dynamics onto a suitable low- dimensional subspace 
(see Section IIV Al for how a low dimensional representation of the macroscopic spatial dis- 
tribution, /io(x, t), is constructed). At each CI step, the coarse- integrated solution is lifted 
from the lower dimensional space into the full microscopic state space and evolved according 
to the Monte Carlo scheme, allowing the error incurred during the coarse time step to relax 
to the slow manifold parameterized by the cell density. Because the internal state {u, s} 
of cells, and their directions of motion, d, are ignored in obtaining fiQ{x,t), they must be 
suitably reinitialized to construct the initial condition for the next MC step (see Section 

nvB}. 
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FIG. 1: Schematic illustration of the coarse integration procedure. 
A. SVD-based restriction scheme 

Below, we outline our coarse integration scheme which operates on data from singular 
value decomposition of the distribution of cell positions from Monte Carlo time evolution 
results. In the following, the solution at a "reporting step" refers to the sorted cell positions 
at specified time intervals, Tgtep- Note that this time interval should be distinguished from 
a single Monte Carlo iteration, of which it contains a large number. 

1. The Monte Carlo description is simulated for k reporting steps, corresponding to a 
total time interval equal to Tgettie- Step by step results are accumulated for a further 
m reporting steps, corresponding to a total time interval equal to Tf^. 

2. At each reporting step, we sort the cell positions, Xi{j), to obtain Aj(j), i = 
k + 1, . . . , k + m, and j = 1, . . . , Aceiis- The matrix Ai is constructed (columnwise): 
Ag = {Xk+i, Xk+2, ■ ■ ■ ,Xk+m}- Figures El and El show the chemoattractant profile, and 
representative time sequences of sorted cell positions and corresponding histograms of 
the spatial density distributions, respectively. However, we find it more convenient to 
work with the cumulative probability distribution function, defined as 

P{X,t) = i /" fioix',t)dx', (17) 

which is related to the sorted cell positions at that time according to 

3lN,Ms = P{X,{j),U), (18) 

where L = Xmax — Xmin is the spatial domain. We will use approximations to P(A, t), 
or equivalently the sorted cell positions, Aj, as the macroscopic variable [26l |. 

3. A low dimensional representation of the sorted cell positions is constructed in terms 
of orthonormal numerical basis functions, {«'•''•'}, 

aM(t,) = «M ■ A, = '5;'«M(j)X,(j), (19) 
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FIG. 2: Chemoattractant profile used in Monte Carlo simulation of chemotaxis: S{x) = So{x ~ 
Xmin)'^{x - Xmaxf, where ^0 = 1.6 X 10~^ and {xmin,Xmax) = (1, H)- The choice of this profile 
ensures zero gradient of chemoattractant at the boundaries, consistent with no-flux boundary 
conditions for the motion of bacteria. 

where the basis functions are obtained from singular value decomposition of the MC 
data as described below. 

We distinguish between two approaches to obtaining the numerical basis set. 

• Local basis functions: The singular value decomposition of A^, constructed in 
Step 2, is computed: 

A, = C/,iy,y/ (20) 

The columns of C/^, given by {ui^^^}, r — 1, . . . , m, are the numerical basis func- 
tions. Wc assume that they remain valid basis functions during the projected 
step. We is a diagonal matrix, where {wei,We2, ■ ■ ■ , w^m} ^'I'c the singular values. 

• Global basis functions: Using a single MC evolution from initial conditions to 
steady state, Tf, the matrix Ag is constructed from the full data according to Step 
2: Ag = {Xi,X2, . . . ,Xm}, where M = Tf/Tgtep- Singular value decomposition 
of Ag gives the the global basis functions, {ug^^^}. (In practice, initial data from 
t — to t — Tgettie are not included.) 

4. We perform linear least squares extrapolation of {a^'^\ti)} using i = k + 1, .... k+m+1 
to obtain {a^'^\tk+m+i+p)} , corresponding to a projected time equal to Tproj- The 
projected solution is given by: 

Yk+i+m+p = ^ a^''^ {tk+m+l+p)u^''^ ■ (21) 
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FIG. 3: Sorted cell positions (top), and corresponding histograms (bottom) using a bin size equal 
to 0.1, at T = (thin solid line), T = 20,000 (dashed line), T = 60,000 (long-dashcd line) and 
T = 100,000 (dot-dashed line), and T = 200,000 (thick solid line). From the time evolution of 
the variance of these distributions we have determined that at T = 200, 000, the coarse solution is 
very close to the stationary state. 
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5. In practice, we use only the numerical basis functions associated with the dominant 
singular values, r = 1, . . . , Vmax, where the truncated number of basis functions, Vmax, 
depends on the spectrum of singular values. 

The use of empirical orthogonal eigenfunctions, also referred to as the Karhunen-Loeve 
expansion or the Proper Orthogonal Decomposition, obtained through singular value de- 
composition of experimental or simulation data, for model reduction in systems with spa- 
tiotemporal dynamics originates with Lorenz in the context of weather prediction 22] and 
has found widespread use in the dynamical systems context since the mid-80's (see the 
monograph and references therein). 

The SVD process is computationally expensive, so if the computational cost of the mi- 
croscopic cell simulations were small, using global basis functions generated in a preliminary 
run would be less costly than using local basis functions. However, in a simulation of a large 
number of cells, the total microscopic simulation cost for the cell population may be much 
larger that for the SVD process. For example, in the numerical experiments reported here 
the CPU time for one simulation step of an individual cell was 2.07 microseconds, while the 
time for one SVD step was 22.4 milliseconds. However, in a simulation of 10^ cells using 
local basis functions, the microscopic integration routine was called 4 x 10^ times versus 16 
calls on the SVD routine so that 0.0041% of the time was spent on SVD. 



B. Reinitialization of internal variables 

The question of consistent schemes for combining microscopic and macroscropic descrip- 
tions of a physical process is central to multiscale modeling. Here, to alternate each CI step 
with MC, in principle allowing relaxation of the numerical error in the projective step, we 
must define an appropriate reinitialization procedure: We know the position of each cell, x, 
after the CI step, but we discard all information about its internal state, {u,s}, given by 
the values of the signal transduction variables and the flagellar states, and its direction of 
motion, dj. 

Our reinitialization protocol is empirically motivated: the signal transduction variables 
are set equal to their local equilibrium values, Uj = {0, f {S{xj))}, where j refers to the 
j*^ cell. This allows the subsequent response of each cell to be within the most sensitive 
range of the motor gain. For simplicity, all (six) flagella are restarted in the CW state, 
Sj = {0, 0, . . . , 0}. In Figure |31 we show relaxation of the ratios of the numbers of right, 
left and tumbling cells at a set location along the chemoattractant profile. These ratios at 
steady state for a flat chemoattractant profile are 

= S, (22) 

where the right and left moving ratios are equal. For the present choice of numerical pa- 
rameters, B ~ 0.86 and Pr = Pl ^ 0.43. Other flagellar reinitialization schemes lead to 
similarly rapid relaxation rates. 
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FIG. 4: Ratios of right moving, left moving and tumbling cells at x = 4 (±0.1), using 
j = 1, . . . , Nceiis = 100, 000 and starting from an initially uniform distribution of cells along the 
chemoattractant profile: Averages over 5 Monte-Carlo realizations started from different random 
number seeds are plotted. The error bars, corresponding to the root-mean-square of the distribu- 
tion of these ratios, are shown for select points. We note fast relaxation of these ratios to their 
equilibrium values. The signal transduction variables were initialized to their local equilibrium 
values, Uj = {0, / {S{xj))}. 

V. NUMERICAL RESULTS 

A. Low-dimensional representation 

In Figures El and IHl we show the first four dominant global basis functions and singular 
values obtained from MC evolution of Nceiis = 10^ from initial conditions to steady state. 
Figure [7| shows the mean coefficients over A^mc = 10 MC realizations 

A''mc 

«^^^W = 7^E4^^W' (23) 
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and the root-mean-square of the distribution of coefficients, given by 




Hence, the error bars in this figure denote the expected statistical variation in the values of 
these coefficients as a function of time. 




FIG. 5: {4 }, = l,---,4 are the first four global SVD basis functions for the cumulative cell 
density. 

In Figure IHl we show the low-dimensional representation of the coarse-integrated solution 
for Tsettie = Tfit = 5 X 10^ and Tproj = 5 x 10^. For each CI step, after evolving the 
Monte Carlo for Tgettie, we use a linear fit to {a^^\t) , a^'^\t) , a^^\t) , a^^\t)} during the 
interval Tfu to project the solution forward in time by Tproj- The points correspond to 
averages over A^mc = 5 realizations, and the error bars give the rms of the distribution 
of these coefficients. For reference, we have included the average MC coefficients in this 
figure. These results indicate that for the higher order coefficients, [a^^\t) , a^'^\t)^ , whose 
coarse dynamics are described by a shorter characteristic time-scale, a higher order time 
integration scheme would be more effective than the explicit Euler method used here (29| . 
Consequently, although for intermediate times the difference between the CI and MC results 
the straightforward CI scheme has captured the macroscopic dynamics of the solution. 
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FIG. 6: Largest singular values from global SVD of Monte Carlo data from initial conditions to 
steady state (Tf = 200, 000, T^tep = 200). The first four dominant modes {vmax = 4) are used to 
construct the low-dimensional representation of the cumulative cell density distribution. 

These results also illustrate the implicit assumption of separation of time scales: The 
coefficients |a^''^(t)}, governing the macroscopic behavior of the system, vary on a time 
scale of (9(10'') units, while the longest microscopic time scale (adaptation time of the signal 
transduction model) is C(IO^). 



B. Analysis of Errors 

For comparison of solutions at different values of the CI parameters, (Tsettie,Tfit,Tproj), 
with the MC, we construct the following measure of relative error: 



Eti M-)(t)-a(-)(t)]^ 



1/2 



(25) 



The coefficients {a^^\t) , a^'^\t) , a^^\t) , a^^\t)} , used in computing the relative error of the 
CI solution are obtained as inner products of the solution with the global basis set, regardless 
(a) of whether global or local basis sets were used in the restricting/lifting step, and (b) the 
dimensionality of the restricted space. 

In Figure ini we first plot the average relative error for A^mc = 10 MC runs, given by: 



MC 



Nmc 

5Z^fc(^)- 



(26) 



k=i 
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FIG. 7: Coefficients of SVD basis functions: a^^'^U) = n^'^) • Xi = Ylf=i"' u^''\j)Xi{j), where the 
points represent averages over A'^mc = 10 MC reahzations, starting from different random number 
seeds. The global basis functions {u^^)} are obtained from one of the MC runs. The error bars, 
plotted for select points, correspond to the root-mean-square of the distribution of coefficients. 



The error bars denote the error on the mean, given by ag = o^j V^mc, where is the 
root-mean-square of the distribution of the error 

f . Nmc ^ 

^s{t) = Yl [^'^W - ■ (27) 

t^MC-l J 

We similarly compute the relative error for CI solutions, obtained using either (i) global 
or (ii) local SVD basis functions in each step. Figure ^1 shows the mean errors for solu- 
tions with CI parameters given by (Tsettie,Tfit,Tproj) = 10^ x (|, 1,5), 10^ x (1,1,10) and 
10^ X (2,2,20). Note that for these parameter values, the efficiency of the CI scheme re- 
mains the same. In this figure, solid lines show results obtained using global SVD basis 

functions, ^^u^g\uf\uf\u^^''^, in the CI steps. Dashed lines show these results obtained 
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FIG. 8: Coefficients of SVD basis functions: cJ'''\ti) — u,g ■ — u,g 

10^ X (i i 5). The points represent averages over A^mc 



CI parameters (T. 



settle 1 ^ fit 1 ^ proj / 



{j)Xi{j), with 
5 

reahzations, starting from different random number seeds. The error bars, plotted for select points, 
correspond to the root-mean-square of the distribution of coefficients. The average MC coefficients 
(dashed line) are shown for comparison. 



using local SVD basis functions. In the latter case, since singular value decomposition is 
applied to a short segment of the MC, we found the higher order local basis functions to be 

relatively "noisy"; hence, were used in each CI step. Note that the rms of the 

distribution of errors obtained for CI solutions using local basis function is larger; in this 
case, the basis functions themselves are also subject to statistical variations. We find that 
the relative errors using global SVD basis functions are generally smaller than those corre- 
sponding to CI solutions obtained using local SVD basis functions. Using either global or 
local basis functions in the CI step, the smallest relative error is achieved with the first two 
sets CI of parameters, corresponding to projected time intervals short enough to resolve the 
macroscopic dynamics of the higher order coefficients of the low-dimensional representation 
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FIG. 9: Relative error given by ek{t) = 

points correspond to the average error, obtained over k = 1, . . . , 10 MC reaUzations starting from 
different random number seeds. The error bars denote the error on the mean, ag. 



of the solution. 

As a benchmark, in Figure [TT] we show the error obtained for (Tsetue,Tfit,Tproj) = 
10^ X (1,1,0), using global SVD functions {ulj'^}, with r = 1,...,4 (dotted line) and 
r = 1,...,8 (solid line). These results show the error incurred in restricting the spatial 
distribution of cell positions to the low dimensional representation given by {a^'^^}, without 
the contribution from the projective step. First, we note that the maximum error obtained 
from the restriction to {a^c^\ . . . , dg^^} and {ctg^'', . . . , af'^} are comparable, indicating that 
using a higher dimensional restriction is not significantly advantageous here. Secondly, this 
error is comparable to the total error in Figure ITUT b). which includes projective integration. 
Given that the total error increases with a longer projective time step, Tproj, as shown in 
Figure Hn^c), these results point out that the "optimal" Tproj, for which the errors due to the 
projective time step and restriction are separately comparable, is achieved in Figure ITUTb). 

Finally, in Figure ^1 we compare CI errors using two different reinitialization schemes 
and find them to be equivalent. 



VI. CONCLUDING REMARKS 

We have demonstrated, through a coarse projective integration algorithm, that short 
bursts of appropriately initialized microscopic simulations can be used to simulate the macro- 
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FIG. 10: Relative error given by e{t) = {Er=i - V Er=i ["^'^n*)]^} > where 

a^^) (t) are the coefficients of the global SVD basis functions for CI solutions corresponding to pa- 
rameters {Tsettle, Tfit, Tproj) = 10^ X (i, i, 5) (bottom), (10^, 10^, 10"^) (middle), and 10^ x (2, 2, 20) 
(top). The CI solution at each step was constructed using the first four global basis functions, shown 
using solid lines, and using the first two local basis functions, shown using dashed lines. The points 
correspond to averages over Nmc = 5 realizations starting from different random number seeds. 
The error bars denote the error on the mean, a^- 



scopic evolution of the cell density distribution in bacterial chemotaxis. The outer level of 
our computational structure was a traditional, continuum Euler integration scheme; the 
time-derivatives of the cell density field it employed, however, were estimated from short 
time evolution intervals of the Monte-Carlo description of this process, and not evaluated 
from a known, macroscopic expression. This approach leads to a general framework for the 
computer-assisted analysis of systems whose dynamics are given at a microscopic/stochastic 
"fine" level, but for which we require averaged, macroscopic information. 

It is interesting to discuss the benefits and shortcomings of such a procedure. If accurate 
closed chcmotactic equations exist, one should use them instead of the two-tier modeling 
we propose here. In addition to the insight gained from exact or approximate analytical 
solutions, computer-aided time-evolution or bifurcation analysis using explicit equations is 
more efficient than kinetic Monte Carlo simulations. However, if such model equations 
are not available, our hybrid computational approach can become more economical than 
long-time direct Monte Carlo simulation. Furthermore, when one is interested in qualita- 
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{Tsettie,Tfit,Tproj) = 10^ X (1,1,0): After each {Tsettie + Tfu) interval, the MC evolution is 
restarted from its low dimensional representation given by |a^^\ . . . , a^"^-* | (dotted line) or 
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, a 
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I (solid line). 



live transitions or bifurcations of the macroscopic behavior, marginally stable or unstable 
stationary states may be difficult to identify through direct microscopic simulations, while 
coarse timestepping holds promise when combined with traditional bifurcation algorithms 

UMM- 

It appears, therefore, that a modeler would ultimately gain in obtaining quantitative 
computational answers efficiently, but perhaps lose in the fundamental understanding of a 
physical process that macroscopic model equations would offer. Therefore a promising re- 
search direction is to use such algorithms to test the validity of explicit closures that assume 
slaving of certain higher order moments of the evolving distribution to lower order ones. In 
the chemotaxis context, when signal transduction and corresponding motor response of the 
cell to an external signal are taken into account as we do here, macroscopic model equa- 
tions derived systematically from the microscopic description do not exist. The assumption 
implicit our the coarse integration algorithm was that the macroscopic description closes at 
the level of the spatial density distribution. In one spatial direction, for example, it would 
be interesting to learn when this assumption breaks down, requiring separate evolution of 
right-going, left-going and tumbling cell density distributions. 

More generally, we believe that the approach illustrated here for coarse time integration 
of a Monte Carlo description of bacterial chemotaxis, is broadly applicable with different 
types of microscopic simulators describing, for example, Brownian, Lattice-Boltzmann, or 
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FIG. 12: Comparison of relative errors for CI solutions (using local SVD basis functions) cor- 
responding to two different flagellar reinitialization schemes. The signal transduction variables 
were initialized to their local equilibrium values, Uj = {0, / {S{xj))}. In one case, all flagella were 
started in the CW state, shown in the thin solid line; in the second case, all flagella were started 

in the CCW state, with those cells to the left (right) of the chemoattractant peak running right 
(left), shown in the thick dashed line. The points correspond to averages over A^mc = 5 realizations 
starting from different random number seeds. The error bars denote the error on the mean, ag. 

molecular dynamics, and leading to emergent macroscopic dynamic behavior. Such a hybrid 
scheme allows efficient simulation of the macroscopic behavior, and may provide insight into 
macroscopic model equations. 
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